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Quantum vortex structures and energy cascades are examined for two dimensional quantum tur- 
bulence (2D QT) at zero temperature. A special unitary evolution algorithm, the quantum lattice 
algorithm (QLA), is employed to simulate the Bose-Einstein condensate (BEC) governed by the 
Gross-Pitaevskii (GP) equation. A parameter regime is uncovered in which, as in 3D QT, there is a 
short Poincare recurrence time. It is demonstrated that such short recurrence times are destroyed 
by stronger nonlinear interaction. The similar loss of Poincare recurrence is also seen in the 3D GP 
equation 1^. Various initial conditions are considered in an attempt to discern if 2D QT exhibits 
inverse cascades are is seen in 2D CT (classical turbulence). In our simulation parameter regimes, 
no dual cascades spectra were observed for 2D QT — unlike that seen in 2D CT. 



I. INTRODUCTION 

A thorough understanding of vortex motion in turbu- 
lence remains an important challenge [2] . Now in classical 
turbulence, the vortex can not be precisely defined and is 
viewed as a continuous property of the fluid. However in 
quantum turbulence the quantum vortex is a topological 
singularity with discrete quantized circulation. The size 
of the quantized (isolated) vortex is characterized by the 
coherent length in superfluid, a measure of the relative 
strength of the kinetic energy. This makes quantum tur- 
bulence, loosely defined as the dynamics of tangled vor- 
tices in 3D or of regularly /irregularly distributed vortex 
points in 2D ^, a suitable prototype for understanding 
classical turbulence [3l |J , particularly in the very long 
wavelength limit in which the discrete quantization of the 
vortices should become unimportant. Additionally, since 
the turbulence driven by a BEC gas expansion displays 
different characteristics from those of a trapped BEC gas 
[5 , there is much still to be understood in the dynamics 
of QT. 

The evolution of quantum turbulence in a 2D BEC gas 
at zero temperature is well described by the mean field 
GP equation for the one particle wave function ?/^: 



iat^/; = -vV + 5lV'lV, 



(1) 



where g is the nonlinear coupling constant for the weak 
interactions among the ground state bosons, in lattice 
units where % = 1 and m — 1/2. Here we omit the chem- 
ical potential term — /iV^ because this term only adds a 
global phase shift e~*^^ to the wave function. One impor- 
tant feature of the GP system is that it is a Hamiltonian 
system with exact conservation of energy. 

If the phase space dynamics of an energy conserving 
system is bounded then there will be a Poincare recur- 
rence of the initial state. For continuous systems, the 
Poincare recurrence period Tp is typically so extremely 
long that for all practical purposes it can be almost con- 
sidered to be infinite. However, for certain discrete maps 
like the Arnold's cat map, Tp can be unexpectedly very 



short. We shall find that there is a class of initial con- 
ditions for the 2D GP equation for which the Poincare 
recurrence time is on the order of (9 [10^] iterations on a 
512^ grid. This Poincare recurrence time is independent 
of the dimensionality of the dynamics - whether 2D or 3D 
[1 . Moreover the Poincare recurrence time scales with 
diffusion ordering: as the (linear) grid is increased from 
L\ to L2, the Poincare recurrence time Tp increases by a 
factor of {L2/ - again independent of the dimension- 
ality of the dynamics. For a short Poincare recurrence, 
the interaction is required to be much smaller than the 
kinetic energy. Under this condition, the BEC gas is 
highly dilute and the recurrence time can be very short. 



A. Superfluid theory 

The simplest theory for a BEC in the zero-temperature 
limit is given by Lagrangian density for a scalar complex 
field 



- VV^^ • VV^- I (V^V) . (2) 



The invariance of ^ with respect to phase rotation, to 
space and to time translations yields the conservation of 
density, momentum and energy, respectively. The result- 
ing set of hydrodynamic equations are: 



+ V • (pv) = 
p{dtv + V • Vv) = -2pV ( g p 



VP 



(3a) 
(3b) 



where p is the probability density and pv the probability 
momentum: 



p = tp^tjj 
pv = — iIj^ViIj). 



(4a) 
(4b) 



Equation (3a) is identified with the continuity equation 



of fluid; while (3b) is recognized as the Euler equation for 



a barotropic inviscid fluid. Notice that the pressure now 
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consists of both local and non-local terms: the normal 
(local) pressure gp^ ^ and the quantum pressure 

which is non-local. To analyze the energy spectra of 
quantum turbulence, the (constant) total energy 



Et = j dx'(^-p\v\^+9p' + 2\V^p\''^ 



(5) 



(in 2 spatial dimensions) is decomposed into three com- 
ponents: 

kinetic energy: Ek ~ 2 J ^^'^ p\^\^ (^^) 
internal energy: Ej = g J dx^ (6b) 
quantum energy: Eq — J I^a/pP- i^^) 

One important feature of quantum turbulence is its com- 
pressibility. Much of the classical turbulence literature 
has focussed on incompressible turbulence but there is 
sound wave emission [6] in a BEG gas. Therefore, to 
compare quantum turbulence with its classical counter- 
part, we further decompose the kinetic energy into its 
compressible and incompressible components. This is 
achieved by introducing the density weighted velocity 
q = y^pv [71 H] and the Helmholtz decomposition of a 
non-singular vector field. The longitudinal component of 
q contributes to the compressible kinetic energy while the 
transverse component contributes to the incompressible 
kinetic energy: 



compressible: Ec = (27r)^ / dk^ |q 



2 IP. |2. 



(7a) 



incompressible: Eic = (27r)^ / dP |qicP; (7b) 
where qc and q^c are defined as: 



- k q 



(8b) 



with q = qc + q^c being the Fourier transform of q. Con- 
sequently, the energy density of incompressible kinetic 
energy and compressible kinetic energy can be defined as 



/•27r 

= k . 
Jo 



dd\q,,Uk,0)\' 



(9) 



using polar coordinates /c, 0. 

For 2D classical incompressible turbulence, the con- 
servation (in the inviscid limit) of both the enstrophy, 
Z = J (ir|Vxi;p, and energy has a profound effect on en- 
ergy transfer. Based on the assumption of incompressibil- 
ity, isotropy and self-similarity, Kraichnan [9 and Batch- 
elor [To] have demonstrated that there exists dual cas- 
cades in the inertial range, with an inverse cascade of 
kinetic energy but a direct cascade of enstrophy in wave 



number space with spectral exponents for the kinetic en- 
ergy spectrum: 



inverse cascade: Sic{k) oc k 
direct cascade: Sic{k) oc k~^ . 



The existence of dual cascades is a hallmark of 2D (clas- 
sical) turbulence. In 3D classical turbulence, enstrophy 
(in the inviscid limit) is no longer conserved and one 
only finds the direct cascade of kinetic energy leading 
to the well known Kolmogorov k~^^^ spectrum. In 2D, 
the inverse cascade dictates that the energy is transferred 
to larger spatial scales, which is manifested by the coa- 
lescence of eddies (with the same rotational sense) into 
larger and larger vortices. On the other hand, the di- 
rect cascade determines that the vorticity is cascaded to 
smaller scales and dissipated away through viscosity at 
the end of the inertial limit. However, for quantum tur- 
bulence, vortices can be created via nucleation [11] and 
annihilated through the mutual interaction leading to the 
fusion of counter-rotating pairs, presumably due to the 
quantum pressure in Eq. (3b). Thus enstrophy in 2D 
QT is not conserved. Lacking the conservation of enstro- 
phy and because of compressibility effects, 2D quantum 
turbulence does not necessitate the existence of dual cas- 
cades. 

The recent papers on 2D quantum turbulence by Horng 
et al. [12 and Numasato et al. [8 , have examined the 
energy cascades. Horng et. a.[12j consider the initial 
conditions of a Gaussian BEG cloud with embedded vor- 
tices in an external strong potential well. In their quasi- 
stationary turbulent state they see a clear demarcation 
between the spatial distribution of the compressible ki- 
netic energy from that of the incompressible energy. In 
particular, the quantum vortices are predominantly lo- 
cated outside the Thomas Fermi radius along with the in- 
compressible kinetic energy distribution, while the com- 
pressible kinetic energy distribution is predominantly lo- 
cated within the Thomas Fermi radius. They find a dual 
cascade of the incompressible kinetic energy spectrum — 
akin to 2D classical turbulence: an inverse energy cascade 
with spectrum /c" with a ~ —5/3 (but with large error 
bars in a very restricted k-range) and a direct enstrophy 
cascade with (a noisy) spectral exponent a ~ —4. This 
deviates from the standard 2D classical direct enstrophy 
cascade exponent a = —3. Numasato [8"^, on the other 
hand, have no external trapping potential and the ini- 
tial state is spatially homogeneous with random phase in 
momentum space. In contrast to Horng et al., they find 
a direct cascade of incompressible kinetic energy with a 
transient spectral exponent a ~ —5/3 (also with large 
error bars and quite a narrow time window). In their 
final asymptotic state, the vortices disappear from their 
system with the incompressible kinetic energy tending to 
zero. 
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B. Organization 

This paper is organized as fohows. In Sec. |TI| we give 
only a brief description of the QLA algorithm since the 
3D algorithm has been given in considerable detail else- 
where [13 . In Sec. |III[ we observe the Poincare recur- 
rence under two very different sets of initial conditions: 
in the first case when there are initially quantum vortices 
present in the system, and in the second case where the 
initial density is uniform but the wave function has ran- 
dom phase fluctuations. Thus in this case, there are ini- 
tially no vortices, but the high kinetic energy will rapidly 
lead to their creation. By gradually increasing the ratio 
of the interacton energy to the kinetic energy, we demon- 
strate the destruction of this fast Poincare recurrence. 
In Sec. [rV] we examine the power law of the energy cas- 
cades. A power law in the incompressible kinetic 
energy is observed, while the compressible kinetic energy 
spectrum is somewhat more complicated and, as in 3D 
QT , exhibits multi-cascade behavior. However, this in- 
compressible power law is absent during the evolution 
of the turbulence when the vortices are annihilated from 
the system and it reappears when vortices are re-created. 
It is suggestive to correlate the behavior of the large k- 
range with this power law spectrum to that found by 
Nore et al. [7] on taking Fourier transform of an isolated 
quantum vortex. Of course, there are some differences: 
for the isolated quantum vortex there is no compressible 
kinetic energy, while from our 2D QT simulations in the 
large k region the compressible kinetic energy spectrum is 
also k~^ and of similar magnitude to the incompressible 
spectrum. 

A close examination of the incompressible energy spec- 
tra, in our chosen parameter regimes, does not reveal the 
existence of the Komolgorov k~^^^ spectrum which char- 
acterizes the inverse energy cascade in 2D classical in- 
compressible turbulence. In Sec. |V| we summarize some 
of our conclusions. 



II. QUANTUM LATTICE ALGORITHM 

We shall find that the GP system can be recovered 
from a simple two qubit QLA algorithm with basis set: 
|00), |01), 1 10), 1 11) at each grid point. Applying a series 
of interleaved local unitary collision and unitary stream- 
ing operations at each lattice site, the QLA algorithm 
models the dynamics of GP system in the long wave- 
length limit on taking moments. To recover the GP for 
the single particle wave function |T3l [14] one need only 
consider just two of these basis set elements. Conse- 
quently, the collision operator C and streaming operator 
S can be reduced to 2 x 2 matrices. The VSWAP gate 
is chosen as the collision operator in our simulation: 



with = I, the identity operator. The streaming oper- 
ator is simply a shift operator defined by: 

5A^.,o = n + e^^-^-*n; (11) 
5a..,i =^ + e^"'^^-n; (12) 

where n = ~ ^ ~ 2^"*" ^ stan- 
dard Pauli spin matrix. Ax^ is the displacement along 
the lattice direction. The subscript a in S/\xi,a indi- 
cates the particular qubit being streamed. Interleaving 
the non-commuting operators C and 5, one obtains the 
unitary evolution operator U for 2D GP system ^ : 

CSAx,,aC. (14) 

The GP wave function is recovered from the qubit rep- 
resentation by il^{v^t) = go(r, t) +(7i(r, t), where qjs is 
the (complex) probability amplitude of |01) and |10), for 
/3 = or 1 respectively. On applying this evolution op- 
erator to this spinor state: 

\q,{t^At))-^-\q,(t))^ 
one obtains the dynamical evolution of the wave function 

m 

i!{v, t + At)= 4>{r, t) + iAa;2 i vV(r, t) + {-lYO[Ax\ 

(16) 

on performing a Taylor expansion with respect to Ax. 
Under diffusion ordering Ax^ ~ At, this equation re- 
covers the free particle Schordinger equation up to order 
0[Ax^]. (It should be noted that the equation itself is 
ostensibly of order (9[Ax^] - so that the numerical error 
term is of order (9 [Ax].) Now the error terms take the 
opposite sign for the evolution operator affiliated with 
the different spinor component. Hence by applying both 
evolution operators on the spinor state q= ( ) : 

q{t^At)=UiUoq{t), (17) 

the error itself can be further reduced to order 0[Ax'^]. 

To incorporate the effect of a potential, we introduce 
the following unitary operator: 

n[V{r,t)]=e-'^'^^''^'\ (18) 

V{r) will be defined to be the nonlinear term in GP equa- 
tion: ^IV^p. With this implementation, the GP equation 
can be reproduced up to the order (9[Ax^] through the 
following unitary operation: 

tlj{t + At) = qo{t + At) + gi(t + At); 

qit + At) = UiQ[V{t + At/2)/2]Uon[V{t)/2]q{t), (19) 

where the potential V{t + At/2) is computed from the 
tlj{t + At/ 2) obtained via the one evolution operation 
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UQQ.\y{t)/2]. Note that diffusion ordering is used in ob- 
taining GP equation from QLA. This adds one additional 
parameter a = At in the GP equation to reflect the spa- 
tial/temporal mesh of the lattice. Consequently, the GP 
equation simulated by Eq.(19): 



(20) 



Due to the unitarity of QLA algorithm, the norm of 
the spinor is automatically conserved. However, in the 
standard QLA algorithm there is a very small numerical 
loss in the mean density during the simulation. This is 
largely due to the overlapping of the two components of 
the spinor: 

5-p = j dv(\^\^-\q\^) = |^r(|go + gip-kop-kiP) 

= j dr{q',q,^qoqt). (21) 

If qoqi is kept purely imaginary during the simulation, 
the overlap between the two components vanishes and 
the mean density is exactly conserved. To achieve this 
goal, two modiflcations are introduced to our previous 
QLA algorithm: 

(1) the initialization of the spinor is chosen to be 
qo qi = 

(2) the potential operator Q is replaced by the non- 
diagonal matrix 



-ia^VAt _ 



COS At] 
-ism[VAt] 



-ism[VAt]\ 
cos[VAt] J ' ^^^^ 



such that > {Qn 



, ^uujv — v3 '^'^^^ip in order to reproduce 

7 

GP equation. 

The QLA algorithm consists of three main compo- 
nents: collision, streaming and local phase rotation to 
introduce the potential. It is clear that streaming does 
not alter the phase of go ^i, e.g. q^qi will remain purely 
imaginary. Since the collision operator C is the v^SWAP 
gate, then applying this collision operation four times 
results in the identity operator: 



(23) 



where a and h are real. In actually applying the full 
QLA algorithm, we note that the collide-stream unitary 
operators do not commute. However, what is critical is 
that the streaming operator never changes the phases of 
the spinor components so that 



cs-^cscs-^cs 



(24) 



of the spinor: 



N 



ib 



acos[VAt]^bsm[VAt] \ 
-iasm[VAt] ^ibcos[VAt]J ' 



As a result, the averaged density 



-iVAt 



{a^ib)f 



(25) 



(26) 



is conserved exactly. 



III. POINCARE RECURRENCE 

To probe the Poincare recurrence in the 2D GP sys- 
tem, we introduce the time averaged ratio of the internal 
energy to the kinetic energy deflned in (|6| as 
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ML 

EK{t) 



(27) 



We shall flnd that 7 plays a crucial role in determining 
the existence of short Poincare recurrence times. Another 
useful quantity is the density weighted vorticity ^ : 



(28) 



with p and v being deflned in Q. The two components 
of uOq are: 

• (^a/p) ^ ' which reflects the variation of the 
density near isolated vortex cores or near entangled 
vortices (i.e., vortices that strongly interact with 
each other); 

• s/pV xv-Gz = y^F ^(r — To), which pinpoints the 
2D location of the vortices with their circulation F. 

Another convenient way to visualize the vortices will be 
to identifying the branch cuts and singularities in a phase 
plot of the wave function ip. In particular, at the location 
of a vortex there will be a phase change of 27r per winding 
number, giving the impression of a branch point and the 
consequent branch cut. To reflect the variation in the 
total number of vortices, the density weighted enstrophy 



-I 



dr \uja 



(29) 



is also calculated in the simulation. 

The color scheme for the distribution plots is "ther- 
mal" : blue stands for the low values while red represent 
the high values. 



Vortices initially imbedded in a Gaussian BEC 
cloud 



for some real e and /. Finally, the modifled potential 
operator Q n does not alter the phase of either component 



For the flrst set of numerical simulations, we assume 
an initial wave function with vortices embedded in an 



inhomogeneous Gaussian BEG background. The total 
angular momentum is chosen to be zero and periodic 
boundary conditions are enforced. To satisfy these two 
constraints, the total wave function is taken as the prod- 
uct of the vortex wave functions: ni'^il^)^ where the 
wave function of each vortex tpi takes the form 

4 

V^tot(r, 0) = h e-^^^ Y[^,{r-r,) (30a) 

i=l 

tlji = tanh(va|r - |)e±^^^^^(^-^^). (30b) 

|'0i(ri)| = at the vortex core itself, n is the winding 
number and we first consider 4 vortices embedded in the 
Gaussian background, h controls the wave function am- 
plitude and a can be viewed as a spatial rescaling pa- 
rameter to resolve flow structures in the turbulence. In 
Fig. |3] and |4] we show the initial conditions for winding; 
number n = 1 and winding number n = 2. The param- 
eters are h 0.05, a 0.01, Wg 0.01, g 5.0. The 
distance between neighboring vortices is: L/4 for n = 1, 
and 2I//11 for n = 2 winding numbers. From the vor- 
ticity and phase plot, the location of the vortices are im- 
mediately pinpointed. The blue dots in enstrophy plots 
indicate clockwise rotating vortices [white circles in phase 
plot Fig.[3](m)], while red dots indicate counter-clockwise 
rotating vortices (black circles in the same phase plot). 
In essence the winding number n = 2 vortex is a doubly 
degenerate simple vortex, with phase change 47r. Grid 
size is 512^. 

The energy ratio parameter 7 = 0.018 for winding 
number n = 1 vortices, and 7 = 0.0036 for winding num- 
ber n = 2 vortices. The total energy is conserved to 7 
significant digits in the simulation for 50000 iteration. 

For such low internal to kinetic energies, we find very 
short Poincare recurrence of the initial conditions. This is 
very evident in the time evolution of the internal energy, 
c.f. Fig.[lJ with the recurrence of the strong initial peak. 
Based on the time evolution of internal energy, it is most 
instructive to examine the spatial distribution of the wave 
function amplitude \iptot\^ the density weighted vorticity 
ujq and the phase at times t = 10500, t = 21000, t = 
31300 and t = 41900, where Tp = 41900, Fig. [3| A 
regular lattice of vortices only occurs at integer multiples 
of the half-Poincare time, mTp/2. For example at t = 
8000, Fig.[2j no such self-similar structure is seen. This is 
characteristic for all times away from integer mTp/2. Of 
course, there is an overall 4 quadrant symmetry because 
of the intial symmetry which is faithfully preserved by 
the QLA algorithm throughout the run. 

At t = 10500 and t = 31300, 16 vortices are present in a 
highly symmetric lattice pattern, which is simply a four- 
fold version of the initial wave function. At t = 21000, 
the distribution of \tl^tot\ and ujq would be same as the 
initial condiiton if the origin of the domain was shifted 
to To ^ To + L/2ex + L/2ey with L being the domain 
size. At t = 41900, the Poincare recurrence period Tp, 
the distribution of both l^^l and ujq closely approximates 
the initial state except for some sound wave interference 
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FIG. 1. The time evolution of the internal energy for winding num- 
ber (a) n = 1 and (b) n = 2 vortices. The first peak corresponds 
to the semi-Poincare period Tp/2 and second peak corresponds to 
the full Poincare period Tp, where Tp = 41900 for a spatial grid of 
512^. The fiuctuations in the internal energy are stronger for n = 2 
vortices, indicative of more vortex annihilations and creations for 
the higher winding number case. 



particularly near the boundaries. In the spatial distribu- 
tion of the phase information there are phase shifts S 
at the branch points/point vortices at Tp. In particular, 
one finds a phase shift of (5 = 7r/4 for counter-clockwise 
rotating vortices while having a phase shift S = — 7r/4 for 
clockwise rotating vortices, thus yielding different branch 
cuts connecting the same branch points. 

Initial vortices with winding number n = 2 are energet- 
ically unstable and will rapidly split into two n = 1 vor- 
tices [16]. Therefore the Poincare recurrence should now 
be viewed as the reproduction of that state immediately 
following the degenerate vortex splitting. To compare 
and contrast the evolution of the initial winding number 
n = 2 vortices with the dynamics of the Poincare recur- 
rence for vortices of winding number n = 1, we again 
consider the dynamics at t = 10500, 21000, 31300, and 
41900. Fig. |4] illustrates the global similarities and lo- 
cal differences with the case of winding number n = 1 
vorticies. Fig. [3] . Around t - 10500 and t - 31300, a 
vortex lattice again is populated by 16 pairs of vortices. 
At the semi-Poincare period Tp/2, the wave function re- 
sembles the initial state but with the origin being shifted 
to (L/2, 1//2). It is assuming to point out that the Arnold 
cap map also exhibits this initial condition mirror sym- 
metry at Tp/2 with full Poincare recurrence at Tp - but 
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(a)lV^totl, t=8000 (b)vorticity uoq, (c)phase 6>, t=8000 
t=8000 




JtW 2jyO 3O0 iOO 5£X) 

(d)line plot near vortex core. 

FIG. 2. The spatial distribution of the amplitude, vorticity and 
phase at t=8000. A clear symmetry is observed between quadrants 
which is due to the periodic boundary conditions and the symmetry 
in the initial conditions. However, within each quadrant, no vortex 
lattice is formed, which is most evident in the amplitude plot (a), 
the vorticity plot (b) and the phase plot (c). Plot (d) is a ID plot 
of l-^l along the y-axis dX x = 197, as indicated by the black dashed 
line in the phase plot (c). This line passes through a vortex pair 
(as identified by l-?/^! =0) . The solid lines in (d) show r near 

the core region. The slopes are: —3.44 x 10~^ for the red line and 
3.46 X 10 "'^ for the green line, as expected from symmetry 



only for certain initial conditions wiki For other ini- 
tial conditions one does not see this mirror symmetry at 
Tp/2. At the full Poincare period Tp, one recovers the 
initial degenerate vortex split state but with consider- 
ably higher background noise than for the n = 1 winding 
number case. This increased noise level can be attributed 
to the higher density of vortices. 

We make a very interesting observation around t = 
10500 for the initial winding number n = 2 vortices: we 
find that a pair of counter-rotating vortices are generated 
between neighboring vortices with the same rotation, c.f. 
Fig. [5j Now in the spatial region between neighboring 
vortices with the same sense of rotation, the flow under- 
goes strong shear. When this shear becomes sufficiently 
strong and exceeds a critical value new vortices are spun 
off by the Quantum Kelvin-Helmholtz (QKH) instability 
[llj. Our simulation confirms the recent postulate [SI [17] 
that the quantum Kelvin-Helmholtz instability can be 
one important mechanism for quantum vortex genera- 
tion. 



B. Random phase initial condition 

We now consider an initial wave function with uniform 
density p but with random phase 6, i/j = ^/pe'^^ . The 
evolution of such a wave function under the GP equation 
is energetically unstable. Randomly distributed vortices 
will be generated rapidly to form turbulence. To gener- 
ate the random phase throughout the lattice domain, a 
bicubic interpolation [18] algorithm is employed. In this 
interpolation the desired 2D function fix^y) is approx- 
imated by the polynomial y) defined on a (normal- 
ized) unit square: 



P{x,y) = ^^ai^jx'y^ 



(31) 



i=0 j=0 



The 16 unknown coefficients ft-j J are determined by en- 
forcing continuity at the 4 corners of the unit square. 
Typically, the continuity conditions are chosen to ensure 
that f{x,y), da:f{x,y), dyf{x,y), da:,yf{x,y) be continu- 
ous at these 4 corners. This yields the required 16 equa- 
tions to determine the coefficients aij uniquely. Thus, 

• Discretize the domain into m x m unit squares, m 
is known as the fragmentation level, with greater 
phase fluctuations for higher m; 

• Generate 4 pseudo-random numbers at each corner 
of the unit square, giving 16 random numbers on 
the unit square; 

• Compute the coefficients aij with the given 16 
random numbers; 

• Periodicity is enforced by equating aij at the 
boundaries of the domain. 

A random phase initial condition, shown in Fig. |6j is 
thereby constructed with m = 8 on a 512 x 512 domain. 

To probe the occurrence of short Poincare recurrence, 
we perform a long-time integration of our QLA represen- 
tation of GP equation to t = 100000 such that the ratio 
of internal energy to kinetic energy is 7 = 0.00287.. In 
Fig.[7|we show the evolution of the mean kinetic Ek and 
quantum Eq energies. (The total energy Et is conserved 
to 10 significant digits throughout the simulation.) In the 
time evolution of the kinetic energy Ek one sees strong 
spikes at t=21000, t=41900, t=63080, t=83980 ... with 
the spike amplitude decreasing. The corresponding phase 
plots at t=21000 and t=41900 are shown in Fig. ^ 



III A 



At t=10500, which corresponds to Tp/4 in Sec. 
no vortex lattice is observed. This can be explained 
by the fact that the vortices are not symmetrically dis- 
tributed initially (in our current case there were no initial 
vortices). Somewhat unexpectedly, the vortices in the 
system are totally absent at t =2 1000, the semi- Poincare 
period Tp/2. A dramatic decrease in the incompress- 
ible kinetic energy Ejc is also observed, c.f. Fig. [Tj^b). 
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(a)|V^tot|,t=0 




o 



(b)t=8000 



(c)t=10500 



(d)t=21000 



(e)t=31300 (f)t=41900=Tp 



(g)cjg,t=0 (h)t=8000 (i)t=10500 (j)t=21000 (k)t=31300 (l)t=41900=rp 




(m)6>,t=0 (n)t=8000 (o)t=10500 (p)t=21000 (q)t=31300 (r)t=41900=Tp 



FIG. 3. Poincare recurrence for vortices with winding number 1. The self-similar structure at t = 10500 and t = 31300 is most notable in 
the vorticity distribution plot, which is caused by the symmetric distribution of the vortices mimicking the initial condition. At t = 41900, 
the amplitude IV^I and vorticity ujq are nearly identical to the initial state. However, the phases of vortices undergo different shifts: 7r/4 
for counter-clockwise rotating vortices and — 7r/4 for clockwise rotating vortices. Spatial grid 512^. 



At t=41900, the Poincare period Tp, the phase distribu- 
tion is globahy restored to its initial state (other than a 
global shift). Again, the vortices disappear from the sys- 
tem with strong dips in the incompressible kinetic energy 
Ejc- The decrease in the Ek spikes can be attributes 
to the loss of fine local structure in the phase 0{t). At 
other instances, we observe random distributed vortices 
upon which the quantum turbulence is developed. Fig. |9] 
illustrates the vortices distribution between t = 20600 
and t = 20800, right before their total disappearance at 
t = 21000. At t = 20600, a large number of randomly 
distributed vortices can be clearly detected which char- 
acterizes the 2D quantum turbulence. At t = 20800 the 
number of vortices is sharply reduced. 



C. Loss of Short Poincare recurrence 

Short Poincare recurrence times occur when the ra- 
tio of internal energy to kinetic energy is very small: 
7 ~ O[10~^]. However, as this ratio increases, the effects 
of the short Poincare recurrence are weakened. With 
the random phase initial condition, our simulations show 
that for 7 ~ O[10~^], the short Poincare recurrence can 
no longer be observed. For t > , many randomly dis- 
tributed quantum vortices are formed. As can be seen 
from Fig. [To) when 7 = 0.0567, there is a strong deple- 



tion of vortices at Tp/2 as the state returns towards its 
initial state of no vortices (see also the strong dip in the 
enstrophy). However, there are no such signatures in the 
energies and enstrophy at Tp. As 7 increases to 0.133, 
no depletion of vortices can be observed, which signals 
the loss of short Poincare recurrence. 



IV. ENERGY SPECTRA OF 2D QUANTUM 
TURBULENCE 

To study the energy spectra, Eq.([9|, in 2D quantum 
turbulence, we first consider the vortex initial conditions 



of Sec. HI A For a 50000-iteration run on a 512^ grid, the 
spectra of incompressible and compressible kinetic energy 
are sampled every 100 iterations. If the energy spectra 
e{k) in Eq§ follow a power law, e{k) oc k then the 
exponent a is immediately retrieved from 9iog(/c) log[£:]. 
Fig. 11 displays the time variation of the incompressible 
spectral exponent Sic for two sets of initial conditions: 
(a) winding number n = 1 vortices, and (b) winding 
number n = 2 vortices. The red horizontal line corre- 
sponds to a k~^ spectrum. To understand the frequent 
appearance of very large spectral exponents for vortices 
with winding number n = 1, we consider the spectrum 
around t ~ 24500. From Fig. 12, it can be readily seen 



that the incompressible energy spectrum undergoes rapid 
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(a)|^tot|,t=0 (b)t=10500 (c)t=21000 (d)t=31300 (e)t=41900 




(k)6>,t=0 



(l)t=10500 



(m)t=21000 



(n)t=31300 



(o)t=41900 



FIG. 4. Poincare recurrence for winding number 2. The wave function evolves in a similar pattern as the n = 1 case. At t ~ 10500 
and t ~ 31300, a lattice of vortices is formed. At t ~ 21000, the wave function restores to its initial distribution with a shifted origin. 
At t = 41900, initial wave function is reproduced. However, more background noise is generated, such as the vortices with small depth 
indicated by the small dots in (h). Spatial grid 512^ 




(a)6>, t=10500 



{h)0, t=10600 



FIG. 5. Quantum Kelvin-Helmholtz instability in a localized re- 
gion on a 512^ grid. The plots depicts the phase distribution at 
t=10500 and 10600 for initial winding number n = 2 vortices. 
The zoomed-in spatial domain is [-256,-128] x [-128,128]. At 
t=10600, a pair of counter-rotating vortices are generated between 
the neighboring vortices with the same rotation, which can be iden- 
tified by the new branch cuts identified by the black arrows. 



change at very small spatial scales (i.e., at large k). This 
irregularity is also demonstrated by the phase plots ^(r) 




FIG. 6. Bicubic fitted initial random phase ^(r, 0) on a 512 x 512 
lattice, with m = 8. The range of the phase is [— 7r,7r], and the 1st 
order derivative of the interpolate is constructed to be continuous. 



at t = 24400, 24500 and 24600 in Fig. 13 One can read- 



ily see that it is the disappearance of vortices which is a 
major cause for the discontinuities in the incompressible 
energy spectrum. 

However, with winding number n = 2 vortices, the ir- 
regular variations in Sic are abated, c.f. Fig. [TT] This 
could be explained by the presence of larger number of 
vortices in the system. With the number of vortices dou- 
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IffE 




(b) 

FIG. 7. Random phase initial conditions, (a): Evolution of the 
total energy Ejp^ the kinetic energy Ex and the quantum energy 
Eq. (b): Evolution of the incompressible kinetic energy Ejc and 
compressible kinetic energy Ec- The spikes in the kinetic energy 
are unmistakably associated with the occurrence of short Poincare 
recurrence. This signature is also observed in (b) with the sharp 
dips in the incompressible energy Ejc- 




(a)phase 6>(r) at (b)t=10500 
t=0 




(c)t=21000 (d)t=41900 



FIG. 8. Poincare recurrence with random phase as initial condi- 
tion. Unlike our earlier case with initial quantum vortices embed- 
ded in a Gaussian cloud, no vortex lattice formation is observed at 
t=10500 (see Fig. [sj. The vortices are randomly distributed and 
resemble nothing like the initial (vortex-free) state. At t=21000, 
the vortices are depleted from the system which now bears some 
similarity with the initial condition. At t=41900, the vortices dis- 
appear again and the global features approximate the initial state. 




(a)vorticity \ujq\, t=20600 




(c)vorticity \ujq\, t=21000 

FIG. 9. The distribution of the (both co- and counter- rotating) 
vortices at times t = 20600 and t = 20800 before their total (but 
momentary) disappearance at t = 21000. At most times in our 
simulation, the distribution of vortices resemble that at t = 20600: 
a large number of vortices randomly distributed throughout the 
whole domain. However, the number of vortices rapidly tends to 
zero as t approaches the semi-Poincare period 21000. 
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(b)Evolution oi Et , Ek , Eq; = 0.133 





(c)Evolution of Enstrophy Z; 7 = 0.0567 



] 2 5 

(d)Evolution of Enstrophy Z; 7 = 0.133 



FIG. 10. Loss of Poincare recurrence as 7 increases, random phase initial condition. The grid size is 512 x 512 for this simulation, 
maximum iteration time is 50,000. The sharp drop of enstrophy in Fig(c) indicates the depletion of vortices from the system, which 
happens around t ~ 21000, the semi-Poincare period. 



Sic 



m im _voo 4eo 500 

(a) winding number n = 1, Sic{t) 



100 w m m 500 

(b)winding number n = 2, Sidt) 



FIG. 11. The slopes of incompressible energy spectra. The fitting range for Sic is /c G [50,100]. The red line indicates Sic = —3. The 
time averaged slopes are: s-ic=-3.23 for n = 1 and Sic = —3.09 for n = 2. The variation of Sic{t) for n = 1 is characterized by irregular 
decrease of Sic, such as the slope near t = 24500. 



bled, it is easier to reach a local critical velocity of counter 
flows to generate new pair of vortices; i.e., QKH occurs 
more frequently. Hence the probability of all the vor- 
tices to be annihilated from the system is extremely low. 
This is not dissimilar to what has been observed in 3D 
quantum turbulence [19 . 

Following [8 , we now examine the effect of random 
phase initial condition to probe the cascades of 2D quan- 
tum turbulence by performing simulations on 32768^ spa- 
tial grids with initial fragmentation level into blocks of 
256 X 256. The total energy is conserved to 6 signifi- 
cant digits for 15000 iterations. The parameters in this 



simulation are: 

coupling constant: g = 1.0; 
initial amplitude: \ip{t = 0)| = 1.0; 

coherence length: ^ = — = 1.0; 

spatial resolution: Ax = 0.03; 
temporal resolution: a = Ax^ = 0.0009; 

size of system: Ls = 32768 x Ax = 983.04^. 

Here we identifying coherence length ^ as is usually done 
in the literature, even though it is strictly defined only 
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FIG. 12. Irregular incompressible kinetic energy spectrum for n = 1. Fitting wave number range: k G [50, 100]. The spectral exponents 
are: (a) Sic = —3.005; (b) Sic = —5.828; (c) Sic = —3.357. The kink (black circle) in (b) suggests sudden changes in the incompressible 
energy at small spatial scales. 




(a)t=24400 (b)t=24500 (c)t=24600 (d)centerleft (e)centerleft (f)centerleft 

blowup, t = 24400 blowup, t = 24500 blowup, t = 24600 



FIG. 13. Phase evolution 0(r,t) for 24400 <t< 24600. The lower panel zooms into a centerleft part of the domain. At t=24500, (e), no 
branch cut can be observed (indicating no vortices present). At t=24400, (d), 6 vortices can be identified, while at t = 24600, (f), there 
are 4 vortices. 



for small perturbations of an isolated vortex fT6^ for a 
specific boundary value problem. 

For a 15000- iteration run, the ratio between internal 
energy and kinetic energy is 7 = 1.55. The time evo- 
lution of the total {Et)^ kinetic quantum {Eq) 
and internal {Ej) energy is shown in Fig.[l4j These time 
evolutions are not dissimilar to that presented by Ref. [8 
with their set of random initial phases. The dynamics 
can be categorized into 2 broad stages besides the usual 
vortex-vortex interactions of quantum turbulence: 

• Rapid generation of vortices. Since the initial wave 
function with random phase is very unstable (because 
of strong velocity field variations arising from v ~ V^), 
a large number of vortices are rapidly generated. This 
generation of vortices causes the originally homogeneous 
amplitude to fiuctuate, which results in rapid 
increase in the internal energy and quantum energy. 
Concurrently, the incompressible energy increases 
rapidly. Provided that the enstrophy also increases 
rapidly at this stage, we conjecture that the increase of 
incompressible energy is mainly due to the generation of 
vortices in the system. 

• Depletion of vortices. In this stage the major energy 
exchange is the energy transfer from vortices to elemen- 
tary excitations, such as long range sound waves. This 
can be confirmed from the decrease of the incompressible 
kinetic energy with a corresponding increase in the com- 



pressible energy with the total kinetic energy remaining 



constant, c.f. Fig 14 



At iteration step t > 3000 (i.e., 0.3 on the figure axis 
with units of 10~^), the only major energy exchange is 
between incompressible and compressible energy. We fo- 
cus on iteration steps t G [4000, 8000] to examine the 
spectra. Fig. [15] depicts the energy spectra at t = 8000 
for both compressible and incompressible energy. The in- 
compressible energy spectra can be broadly categorized 
into four regions: (I) k < 0.01%; (II) 0.01% < A: < 0.1%; 



(III) 0.1% 



^ /c ^ %; (IV) % ^ /c. The separation be- 
tween (I) and (II) is signaled by the sudden drop of com- 
pressible energy at /c ~ 15 (the small dip encircled in 
Fig. 15). The regression fit of the incompressible energy 



spectrum is illustrated in Fig.[T6j The variation of energy 
spectra is continuous in different regions and no bottle- 
neck effect are observed. It is interesting to notice that 
around /c ~ %, the incompressible energy slope a is close 
to —4.0 and persists for a time interval 3000 < t < 14000 
during which the randomly distributed vortices dissipate 
away. 

To further examine the incompressible energy spectra 
near /c ~ %, we sample the incompressible energy spec- 
tra every 50 iteration steps between 6000 < t < 10000 
within the wave number window k G [800,1200]. The 
time averaged slope {a) = —4.145 ± 0.066. This is in 
good agreement with Saffman's k~^ power law observed 
in Horng et. al. [TJ- However, we do not observe any 
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(a)Evolution of Ek, Eq and Ei (h)Ej v.s. Ec 



FIG. 14. Evolution of energies for random initial condition. (a): Evolution of kinetic energy Ex, quantum energy Eq, internal energy Ej 
and total energy Et- (b): Exchange between incompressible kinetic energy (red line) and compressible kinetic energy (blue line). Grid 
32768^ 




10^ r 



1 10 100 lutits ir 

FIG. 15. Spectra for incompressible (red) and compressible (blue) 

27T 

kinetic energy at f = 8000. The unit of momentum k is ku = — , 

Ls 

with Ls = 983.04^. The black dashed line indicates the position of 

27r . . . 

k^ = — . The circle emphasizes the dip in the compressible energy 

which seems to propagate with time to smaller k like a backward 
propagating pulse. Grid 32768^ 

inverse energy cascade for k < k^. 

V. CONCLUSION 

We have investigated 2D quantum turbulence using a 
novel unitary QLA algorithm. The unitarity of the quan- 
tum algorithm makes it an efficient method for simulating 
the dynamics of CP system which demands the conser- 
vation of number of particles and total energy. The local 



collision and streaming operation enables the QLA to be 
parallelized almost ideally. The superlinear paralleliza- 
tion of the QLA allows us to probe the energy cascades 
with large grid simulation, such as 32768^. 

The spectra analysis of the incompressible kinetic en- 
ergy reveals a ubiquitous k~^ power law for large k > k^. 
This power law, in 2D QT, may possibly be attributed to 
the result of Fourier Transform of the topological singu- 
larities (and so may possibly be used to identify vortices 
in the 2D CP system). However, there is still substantial 
sound waves at these very small scales due to the presence 
of a compressible kinetic energy spectrum which overlays 
the incompressible spectrum for k > k^. At /c ~ /c^, a k~'^ 
power law is observed, not dissimilar to that found by 
Horng et. al. [12] who then connect this to the Saffman 
k~^ spectrum in CT. However, no inverse energy cas- 
cade is observed in our simulation. This may attribute 
to the compressibility of CP system and fluctuation of 
enstrophy. During the simulations, we discover an unex- 
pected short Poincare recurrence provided that the ratio 
of internal energy to kinetic energy 7 ^ O[10~^]. When 
7 increases to the order of (9[10~^], this short Poincare 
recurrence is no longer observed. 
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FIG. 16. Regression fits to the incompressible kinetic energy 
spectrum k~^^ for various wave number windows, (a) Region /: 
a = +2.34 for 1 < < 15 (red), Region II:a = +0.65 for 
20 <k < 100 fgreen); (b) Region //: a = +0.65 for 50 < /c < 100 
("green), Region Ill.a = -4.17 for 600 < /c < 1000 (red) ; 
(c) Region ///: a = -4.23 for 700 < k < 1200 (red), Region 
IV.a = -3.03 for 4000 <k < 5000 ("green). Grid 32768^. 



